A general approach to the sign problem 
— the factorization method with multiple observables 
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The sign problem is a notorious problem, which occurs in Monte Carlo simulations of a system with 
the partition function whose integrand is not real positive. The basic idea of the factorization method 
applied on such a system is to control some observables in order to determine and sample efficiently 
the region of configuration space which gives important contribution to the partition function. We 
argue that it is crucial to choose appropriately the set of the observables to be controlled in order for 
the method to work successfully in a general system. This is demonstrated by an explicit example, 
in which it turns out to be necessary to control more than one observable. Extrapolation to large 
system size is possible due to the nice scaling properties of the factorized functions, and known 
results obtained by an analytic method are shown to be consistently reproduced. 

PACS numbers: 05.10.Ln, 02.70.Tt, 11. 15. Ha 



I. INTRODUCTION 

Monte Carlo simulation is a powerful tool for studying 
statistical systems from first principles. When the parti- 
tion function has an integrand which is not real positive, 
however, one encounters a notorious technical problem 
called the sign problem. There have been many pro- 
posals, but most of them are successful only for a very 
special system [l[ or for a very small region of the param- 
eter space [|| . In Ref. Q , two of the authors proposed a 
method termed the factorization method, which is based 
on the factorization property of the density of states of 
properly chosen physical observables. It has been tested 
in a random matrix theory for finite density QCD 4 1 , and 
applied also to the lattice QCD at finite density [f| 6| with 
some important new ideas [6|. In the lattice gauge the- 
ory, a lot of efforts are actively pursued also using other 
methods such as analytic continuation Q , multiparame- 
ter reweighting [8|], complex Langevin dynamics [9( and 
the Taylor expansion method [l(| . 

The basic idea of the factorization method is to con- 
trol some observables in order to determine and sample 
efficiently the region of configuration space which gives 
important contribution to the partition function. While 
the previous studies suggest its potential usefulness 
in a wider range of applications, the choice of the ob- 
servables to be controlled seemed rather arbitrary. In 
this work we argue that it is actually crucial to choose 
them appropriately in order for the method to work suc- 
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cessfully in a general system. With this new insight, we 
consider that the factorization method has become a very 
promising approach applicable to any interesting system 
that suffers from the sign problem. 



II. SIGN PROBLEM 

Let us consider a general system defined by the parti- 
tion function 



dAe -S [A] + iT[A] 



(1) 



where A represents the dynamical variables, and Sq and 
r represent the real part and the imaginary part of the 
action, respectively. Since the integrand of ([T]) is not 
real positive due to T, one cannot view it as a sampling 
probability in a Monte Carlo simulation. One way to 
calculate the expectation value of an observable O is to 
use the reweighting formula 



(0) = 



(Oe' r ) 
<e ir >o 



(2) 



where the expectation values on the right-hand side are 
taken with respect to the phase-quenched model 



Z Q = JdAe~ s °W 



(3) 



which can be simulated in the usual manner. The expec- 
tation value (e ir )o is nothing but the ratio of the partition 
functions Z/Zq. Therefore it decreases exponentially for 
large system size V as e~ VA f , where Af > represents 
the difference in the free energy density of the two sys- 
tems. This can happen due to huge cancellations from 
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e zr , and occurs also in the numerator of (|2]). As a result 
one needs O(e const y ) configurations to compute the ex- 
pectation value (O) with given accuracy. This is called 
the sign problem. 

There is also a general problem called the "overlap 
problem" in using a reweighting formula like @. This 
occurs since the region Rq of configuration space one 
can sample effectively by simulating the phase-quenched 
model has very little overlap with the region R which 
gives important contribution in the evaluation of (Oe* r )o 
and (e ir )o- One has to run a very long simulation until 
one can sample enough configurations in the region R. 



III. FACTORIZATION METHOD 

In order to reduce the overlap problem, we control 
some observables so that one can sample configurations 
in the region R efficiently. Let us introduce a set of such 
observables 



E = {O fc |fc = l,-. - ,n} 
and define the normalized observables 

k (Ofc)o ' 
Their expectation values can be written as 



(0 S ) 



Xj p{X\ 5 * ' * 5 X n ) 



(4) 



(5) 



(G) 



using the density of states 

In \ 
p(x lr -- ) x n ) = ( Y\_S(x k -O k ) 



(7) 



Vfc=l 



Applying the reweighting formula to (J7J, one can easily 
derive the factorization property 

p(xi,--- ,x n ) = ^ p (0) (a;i, ■ • • ,x n )w(x!,--- ,x n ) , (8) 

where C — (e ir )o and 

P (0) (si, - ,x n ) = /f[S(x k -d k )\ (9) 
\fe=i / o 

is the density of states for the phase-quenched model. 
The correction factor w(xi, • • • , x n ) is given by 

w(xi, ■■■ ,x n ) = (e lT ) Xu ... !Xn (10) 

as an expectation value in a constrained system 



~ n 

Z{xi,--- ,x n )= / dAe~ s ° l[S(x k -d 
J k=i 



(11) 



In what follows we assume that w(x±, ■ ■ ■ , x n ) is real and 
positive 11 1 . 

When the system size V goes to oo, the expectation 
values are given by (Ok) = x~k, where (xi, • • • ,x n ) de- 
notes the position of the peak of p(x\, ■ ■ ■ ,x n ). Hence 
they can be obtained by solving 



,%n) = -q ^(^l,' 

OXk 



,X n ) (12) 



for k = 1, 



, n, where we have defined 
1 d 



lim 

y^oo V OXk 

1 



\ogp {0 \x 



,x n ) = lim — logw(xi, ■ 



(13) 



Note that the right-hand side of (Tl2|) . which represents 
the effect of T, can largely shift the peak from that of 
(xi , ■ ■ ■ ,x n ), which is given by x k — 1 due to the 
chosen normalization ([S]). Thus we can obtain (Ok) in- 
cluding the contribution from configurations that are dif- 
ficult to sample by simulating the phase-quenched model 
without any constraints. 



IV. CHOICE OF OBSERVABLES 

What is written above is mathematically correct for 
an arbitrary choice of the set S of observables. We will 
argue, however, that the overlap problem can still occur 
in the evaluation of (flQ)) by simulating (flTT) . Note that 
the only difference from evaluating the denominator of 
(T2|) is the existence of the constraints. In fact it turns 
out to be important to choose the set £ appropriately so 
that the overlap problem is removed. 

In order to clarify the overlap problem in the evalu- 
ation of (TTU|) . let us consider another observable O ra +i, 
and define the corresponding functions p^ and w with 
n replaced by n + 1 . Then we obtain the relation 



w{x\ 



J dx n+1 p (0) (xi, ■ • ■ ,x n+1 ) w(xx, ■ ■ ■ ,x n+1 ) 
J dx n+1 p(°)(xi, ■ ■ • ,x n+1 ) 



(14) 



When one simulates (jlip . one mostly samples configu- 
rations with O n+ i close to (O n +i) xlt ... Xn . However, it 
can happen that the integration over x n+ \ in the numer- 
ator of (|14p has important contribution from the region 

of x n+ i not close to (O n +i)n,-,i„ if w{xi,--- ,x n +i) 
has strong dependence on x n+ i. In that case, one 
clearly has some remaining overlap problem in obtaining 
w(xi, • • • , x n ) by simulating (fTl]) . The overlap problem 
can be reduced further by including the observable O n +i 
in the set £. One can, in principle, repeat this procedure 
until the overlap problem is totally removed. 

Let us then discuss what is the minimal set of observ- 
ables that is sufficient to remove the overlap problem. 
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For that, let E be a set of all possible observables, and 
assume that there is no more overlap problem. Here we 
consider a simplified situation, in which the right-hand 
side of (fT2"]l is nonzero for 1 < k < K, and vanishes for 
the rest; i.e., 
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, x n ) =0 for K < j < n 



(15) 



near the peak (xi, • • • ,x n ). Then we can actually use a 
smaller set E' = {O^ \ k = 1, • • • , K} without having the 
overlap problem in evaluating w(x\, • • • , xk) around the 
peak. Furthermore, from (fT5|) one can easily show that 
the expectation value (Oj) for j > K can be evaluated 

as {Oj) ~ (Oj)^,-.. .x K - That this happens in spite of the 
sign problem can be understood by noting that 



(<^e ir } Sl , 



(16) 



namely the observables Oj (j > K) are decorrelated with 
c' r . The fluctuation of the phase plays an important role 
in the determination of (x±, ■ ■ ■ , xk) through the saddle- 
point equation (|12p . but once they are determined, the 
phase can be neglected completely when evaluating the 
expectation values of the observables Oj (j > K). 

In general, we can expect (|15[) to hold only approxi- 
mately. In that case the systematic error involved in the 
above evaluation of (Oj) is given by (j, I > K) 



Ax-j = (H- 1 ) J i — <S>(x 1 ,---,x ri 
oxi 



(17) 
(18) 



This may also cause a small overlap problem in evaluating 
w[x\, ■ • ■ ,xk) around the peak. 

A practical way to search for the observables to be 
included in the set E is to calculate an ensemble average 
of e* r in the phase-quenched model with an observable 
O fixed to x. If it becomes much larger in some region 
of x, the observable should be considered as a candidate. 
We expect that there are many systems in which only a 
few observables have to be included in the set E. 



V. AN EXPLICIT EXAMPLE 

Let us demonstrate how the method works in an ex- 
plicit example. Here we study a matrix model defined by 
the partition function .12] 

Z= J dAe~ Sb {detV) Nt , (19) 

1 4 4 
S b = -iV^tr^) 2 , V = J2^^^ > (20) 



^=1 



where (// = 1, - • • , 4) are N x N Hermitian matrices 
and the 2x2 matrices are Pauli matrices Tj = Oj 



for j = 1, 2, 3, and — il. The system has a rotational 
SO(4) symmetry corresponding to n- O^A^ with 
O G SO (4). The determinant det V is complex in general. 
Under parity transformation A4 — > —A4, Aj — > Aj (j ^ 
4), it transforms as det 2? —> (detP)*. This implies that 
det£> is real for configurations with A4 — 0. From this 
fact alone, it follows that the phase of the determinant 
becomes stationary for configurations with A4 = A3 = 
since one cannot have a phase fluctuation within a linear 
perturbation around such configurations [i"3j |. 

This model was proposed 12'] as a toy model for the 
spontaneous symmetry breaking (SSB) of the SO(10) 
rotational symmetry expected to occur in the IKKT 
model, a conjectured nonperturbative formulation of su- 
perstring theory [3]. In the IKKT model, the space- 
time is represented by the eigenvalue distribution of A^ 
(M = 1, ■ • • , 10) [IE El, and the SSB of SO(10) realizes a 
scenario for dynamical compactification of extra dimen- 
sions in superstring theory. This scenario is supported 
by explicit calculations based on the Gaussian expansion 
method (GEM) [H-il|. It is considered [13 that the SSB 
is induced by the phase of the complex Pfaffian obtained 
by integrating out the fermionic variables. 

In the model (fl9)) . the SO (4) rotational symmetry is 
expected to be spontaneously broken in the large- -/V limit 
with r = N[/N fixed, where Nf is the exponent in (|19p . 
As an order parameter, we consider the "moment of in- 
ertia tensor" 



Tfiv — — tr (A^A U ) 



(21) 



and its real positive eigenvalues Afc (fc = 1, • • • ,4) ordered 
as Ai > A2 > A3 > A4. If their expectation values turn 
out to be unequal in the large- N limit, it implies the SSB 
of the SO (4) symmetry. 

The model CE]) was studied by the GEM up to the 
ninth order for r < 2 [2^]. It was found that the SO (4) 
symmetry is spontaneously broken down to SO(2). Fur- 
thermore, by controlling the eigenvalues to have a hier- 
archy A<j 3> Xd+i, one can obtain d-dimensional config- 
urations (d = 1,2,3), for which the phase fluctuations 
become milder according the arguments below Eq. (|20|) . 
These properties of the model makes it an ideal testing 
ground for the ideas in Sees. IIIII and HVl 

At r — 1, for instance, the results obtained by GEM 
are (Ai) = (A 2 ) ~ 2.1, (A 3 ) ~ 1.0, (A 4 ) ~ 0.8, whereas 
(A/c)o = f for the phase-quenched model. Therefore, the 
expectation values of observables normalized as ([5]) are 

(Ai) = (A 2 ) ~ 1.4 , <A 3 ) ~ 0.7 , <A 4 ) ~ 0.5 . (22) 



VI. MONTE CARLO SIMULATION 



We study the model (fl9j) at r — 1, where the sign 
problem is severe. Since we know that the eigenvalues 
Afc have strong correlation with the fluctuation of T, we 
use Ok = Afc (k = 1, • • ■ , 4) as the observables in the set 
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E. Searching for the solution to (| 12[) in its full general- 
ity is a formidable task. Here we simply check that the 
GEM result (|22|) is indeed a solution to the saddle-point 
equation ([T2"j). This is sufficient for demonstrating that 
the method is error free and that the generalization of 
the factorization method proposed in this work is neces- 
sary in order to remove the remaining overlap problem. 
The question whether one can successfully determine the 
absolute maximum of p(x±, . . . , X4) is subject to a more 
detailed investigation. 

Let us therefore assume that the SO (2) symmetry re- 
mains, which implies x\ = x%- Equation (|12[) then re- 
duces to 



d 
dx k 



\ogPso {2) {x 2 ,x 3 ,x 4 ) 



d 

dx k 



log W SO ( 2 ) {X 2 ,X 3 , X 4 ) 

(23) 



for k = 2,3, 4, where we have defined 



Pso(2) ( x > z ) = p { "' ^ y> z ) 

wso(2) {%, V, z) = w(x, x, y, z) . 



(24) 



CO 

_o 

■><" 



N=16 — a—i 
N=32 o 
N=64 




FIG. 2: The function ^ ■£ log p (0) (a;) is plotted for N = 
16, 32, 64. The dashed line is drawn to guide the eye. We also 
plot — 4- \imN^aa{-^i log w(x)} obtained from Fig. [T] where 
the two solid lines show the margin of error. The value of x 
at the intersecting point is consistent with the GEM result 
<A 2 ) ~ 1.4. 



In order to be brief [23j, let us only discuss Eq. (|23|) 
for k — 2. We set X3 — 0.7 and X4 = 0.5, and solve the 
equation for x%. We also define w(x) = wso(2)( a; ) 0.7, 0.5) 

and p(°\x) = Pso (2) (x, 0.7, 0.5). 
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FIG. 1: The function log«;(a;) is plotted for N = 
The solid line represents the asymptotic behavior 
the coefficients extrapolated to N = 00. 



In Fig. Q] we plot our Monte Carlo data for log w(x). 
It approaches zero at large x, where the dominant con- 
figurations have A2 3> A3, and the phase T becomes sta- 
tionary due to the argument below Eq. (f20l) . One can 
actually show that the asymptotic behavior of w(x) at 
large x is 



1 

N 



2 \ogw(x) 



-C\X 



C 2 X 



-5/2 



(25) 



Fitting our data to J25J) for N = 8, 12, 16, we find that 
finite N effects in the coefficients c\ and c 2 are consistent 
with 0(l/iV). Making extrapolations to N — 00 based 
on this observation, we obtain c\ = 0.322(2) and c 2 — 
0.021(1). The solid line in Fig. [T] represents Eq. ([2g|) 
with these extrapolated values. Figure [5] shows that the 
solution is x = 1.373(2), which is consistent with the 
GEM result (A 2 ) ~ 1.4. 

Let us then consider what happens if we add 



(26) 



as the fifth observable in the set E. We define the corre- 
sponding functions p^ and w with five arguments, and 
also define the reduced functions 



4 0) < 

wo{x) 



p b } (x) = (1.4, 1.4, 0.7,0.5,3!) , 
u;(1.4,1.4, 0.7, 0.5, x) , 



(27) 



which correspond to fixing Afc to the GEM result (|22j) . 
and O = 0/(0)o to some value x. We find that 
-4? log wo (x) approaches zero as x — > [24j with the 
asymptotic behavior 



— log w (a;) = -d lX 2 + d 2 x 5/2 



(28) 



Therefore, the observable O may, in principle, be a "dan- 
gerous" one which must be included in the set E. Figure 
[3] is a plot obtained similarly to Fig. [3J which shows that 
the effect of the phase represented by ([IT]) is to shift the 
estimate of (O) by Ax = 0.07(3). On the other hand, 
the standard deviation of the distribution p^\x) is es- 
timated as a ~ 0.7 /N from the slope of the function 
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plotted in Fig. [3] around x ~ 0.92. This means that the 
deviation Ax is < 2 tr for iV < 16. Thus, the remaining 
overlap problem in obtaining the data in Fig. Q] is prac- 
tically small as far as this observable (|2"6"|) is concerned. 
The fact that we are able to reproduce the GEM result 
with only four observables Xk in the set £ suggests that 
the remaining overlap problem is indeed not so severe. 



o 



O 



X 



3 
2.5 

2 
1.5 

1 

0.5 


-0.5 
-1 



N=8 — a- 
N=16 <—*■- 
N=32 



(d/dx) (1/N") log w (x) 



0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 

x 



FIG. 3: The function ^ ^ log p$ (x) is plotted for N = 
8, 16, 32. The dashed line is drawn to guide the eye. We also 
plot — -4- limjv-»oo{-j^i2- log wo(x)}, where the two solid lines 
show the margin of error. From this figure, we find that the 
position of the peak shifts from x = 0.92 to x = 0.85(3) due 
to the effect of the phase. 



VII. SUMMARY AND DISCUSSION 

In this work we have discussed a general approach to 
systems with the sign problem based on the factorization 



property |8]) of the density of states. The method aims at 
reducing the overlap problem by controlling the observ- 
ables which have strong correlation with the phase e lT . 
Once the solution to the saddle-point equation (fT2"j) is ob- 
tained for such observables, all the other observables can 
be investigated without the sign problem. We have pre- 
sented an explicit example, in which Monte Carlo data 
suggests that one has to control only a few observables 
to remove the overlap problem almost completely. We 
speculate that this is the case in many interesting sys- 
tems. Finding the minimal set of observables for each 
system as described in Section HVl would be the subject 
of a future investigation. 

The main task of the method is to calculate the func- 
tion pop by simulating (fTT|) . which still suffers from can- 
cellations due to e* r . However, near the solution of the 
saddle-point equation (|T2l) , the fluctuation of T becomes 
milder than in the phase-quenched system without con- 
straints. Hence it is expected in many cases that one can 
compute the function (|10[) . directly or by using asymp- 
totic behaviors such as ([25]) . for reasonable system size. 
Then the scaling properties represented by (fT3|) enable 
extrapolations to infinite system size. We hope that 
Monte Carlo studies of many interesting systems that 
are hindered by the sign problem can be made possible 
by using the factorization method. 
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